In [25]:
%pylab
%matplotlib inline
import astropy.io.fits as fits
from photutils import daofind
from photutils import CircularAperture,SkyCircularAperture
from astropy.stats import sigma_clipped_stats
from photutils import aperture_photometry
import photutils.utils.wcs_helpers
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
from astropy.coordinates import SkyCoord
from astropy import wcs
import astropy.units as u
rcParams['image.cmap'] = 'viridis'
Upload the flattened files to Astrometry.net, download "new-image.fits" with sky coordinates added.
Photutils basics via:
In [26]:
#load example astrometry.net frame:
output_string = "/mnt/camp-storage/ATC2016/transition_disk_accretion_rate/c7560t"
nims = 22
first_frame = 3
basestring = "/mnt/camp-storage/ATC2016/wiyn_09/Converted2016june20/c7560t"
output_file = output_string+'%04i' % first_frame +'flattened.fits'
hdu = fits.open("/mnt/camp-storage/ATC2016/transition_disk_accretion_rate/new-image(1).fits")[0]
#(this should be replaced with a public dataset)
#find bright sources
image = hdu.data
bg=np.median(image[:200,:200])
data=image-bg
mean, median, std = sigma_clipped_stats(data, sigma=3.0, iters=5)
print(mean, median, std)
sources = daofind(data-median, fwhm=10, threshold=7.*std) #subtracting the
#so much noise it takes an oversized gaussian to knock it down and not just get hot pixels
#make sure the sources all look like stars:
plt.figure(figsize=[8,8])
positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=20)
#these apertures are just for checking we found stars
norm = ImageNormalize(stretch=SqrtStretch())
plt.imshow(log10(image), cmap='plasma', origin='lower')#, norm=norm)
apertures.plot(color='black', lw=1., alpha=0.5)
In [27]:
#look at an aperture up close
plt.imshow(log10(image), cmap='plasma', origin='lower',interpolation='nearest')#, norm=norm)
apertures.plot(color='black', lw=1., alpha=0.5)
plt.xlim([1000,1100])
plt.ylim([1000,1100])
Out[27]:
In [ ]:
In [ ]:
In [28]:
#convert the xy coordinates of the sources to RA and DEC
w = wcs.WCS(hdu.header)
source_skycoods = photutils.utils.wcs_helpers.pixel_to_icrs_coords(sources["xcentroid"], sources["ycentroid"], w)
c = SkyCoord(source_skycoods[0],source_skycoods[1],frame='icrs')
#define the size of the photometry apertures
apertures = SkyCircularAperture(c, r=7. * u.arcsec)
In [29]:
imshow(hdu.data,vmax=200)
colorbar()
Out[29]:
In [30]:
#do photometry on each file
import glob
astrometry_files=glob.glob("/home/edouglas/camp-storage/ATC2016/transition_disk_accretion_rate/astrometry/0*.fits")
phot=np.empty([len(apertures.positions.ra),len(astrometry_files)])
for i,file in enumerate(astrometry_files):
print(file)
hdu = fits.open(file)[0]
phot_table = aperture_photometry(hdu, apertures)
phot[:,i]=np.array(phot_table['aperture_sum'])
#this is used above to output the file
In [ ]:
In [31]:
#calculate and plot the difference in magnitudes versus the brightest star
subplot(121)
brightest_row=np.where(phot==phot.max())[0]
delta_mags=(-2.5*np.log10(phot)+2.5*log10(phot[brightest_row])).T
plt.plot(delta_mags)
plt.ylabel("$\Delta$mag")
plt.xlabel("frame number")
subplot(122)
#delete dimmest star
phot_pared=np.delete(phot,where(delta_mags==np.max(delta_mags))[1],axis=0)
delta_mags=(-2.5*np.log10(phot_pared)+2.5*log10(phot_pared[np.where(phot_pared==phot_pared.max())[0]])).T
plt.plot(delta_mags)
plt.ylabel("$\Delta$mag")
plt.xlabel("frame number")
Out[31]:
In [32]:
#calculate the delta mag of the target star versus the sum of the others
diff_phot = -2.5*np.log10(np.sum(phot_pared,axis=0)) + (2.5*log10(phot_pared[brightest_row]))
#this is very simplified, another approach: http://adsabs.harvard.edu/abs/1992PASP..104..435H
plt.plot(diff_phot.T)
plt.ylabel("$\Delta$mag")
plt.xlabel("frame number")
Out[32]:
In [ ]:
In [ ]: